/*
** --------------------------------------------------------------------------
**
** auxFunc.h  -  Auxiliar functions to FE process.
**
** Author - Antonio Miranda, January - 2008
*/


#include <stdio.h>
#include <stdlib.h>
#include <string.h>


/* ============================== CMatrix ==================================*/
void CMatrix (double E1, double E2, double E3, double Nu12, double Nu21, 
              double Nu13, double Nu31, double Nu23, double Nu32, double **C)
{
 double G1, G2, G3;
 int    i;

 for (i = 0; i < 6; i++)
   memset (C[i], 0, 6 * sizeof (double));
 
 G1 =  ( E1 * E2 ) / ( E2 + E1 * ( 1.0 + 2.0 * Nu12 ) );
 G2 =  ( E1 * E3 ) / ( E3 + E1 * ( 1.0 + 2.0 * Nu13 ) );
 G3 =  ( E2 * E3 ) / ( E3 + E2 * ( 1.0 + 2.0 * Nu23 ) );
 
 C[0][0] =  E1 * (Nu23*Nu32-1.0)  / (Nu12 * (Nu21+Nu23*Nu31) + Nu13 * (Nu21*Nu32+Nu31) + Nu23*Nu32 - 1.0);
 C[1][0] = -E2 * (Nu21+Nu23*Nu31) / (Nu12 * (Nu21+Nu23*Nu31) + Nu13 * (Nu21*Nu32+Nu31) + Nu23*Nu32 - 1.0);
 C[2][0] = -E3 * (Nu21*Nu32+Nu31) / (Nu12 * (Nu21+Nu23*Nu31) + Nu13 * (Nu21*Nu32+Nu31) + Nu23*Nu32 - 1.0);
 C[0][1] = -E1 * (Nu12+Nu13*Nu32) / (Nu12 * (Nu21+Nu23*Nu31) + Nu13 * (Nu21*Nu32+Nu31) + Nu23*Nu32 - 1.0);
 C[1][1] =  E2 * (Nu13*Nu31-1.0)  / (Nu12 * (Nu21+Nu23*Nu31) + Nu13 * (Nu21*Nu32+Nu31) + Nu23*Nu32 - 1.0);
 C[2][1] = -E3 * (Nu12*Nu31+Nu32) / (Nu12 * (Nu21+Nu23*Nu31) + Nu13 * (Nu21*Nu32+Nu31) + Nu23*Nu32 - 1.0);
 C[0][2] = -E1 * (Nu12*Nu23+Nu13) / (Nu12 * (Nu21+Nu23*Nu31) + Nu13 * (Nu21*Nu32+Nu31) + Nu23*Nu32 - 1.0);
 C[1][2] = -E2 * (Nu13*Nu21+Nu23) / (Nu12 * (Nu21+Nu23*Nu31) + Nu13 * (Nu21*Nu32+Nu31) + Nu23*Nu32 - 1.0);
 C[2][2] =  E3 * (Nu12*Nu21-1.0)  / (Nu12 * (Nu21+Nu23*Nu31) + Nu13 * (Nu21*Nu32+Nu31) + Nu23*Nu32 - 1.0);
 C[3][3] =  G1;
 C[4][4] =  G2;
 C[5][5] =  G3;
}

/* ============================== InvCMatrix ==================================*/
void InvCMatrix (double E1, double E2, double E3, double Nu12, double Nu21, 
                 double Nu13, double Nu31, double Nu23, double Nu32, double **D)
{
 double G12 = (E1*E2)/(E2*(1.0 + Nu21) + E1*(1.0 + Nu12));
 double G13 = (E1*E3)/(E3*(1.0 + Nu31) + E1*(1.0 + Nu13));
 double G23 = (E2*E3)/(E3*(1.0 + Nu32) + E2*(1.0 + Nu23));
 int i;

 for (i = 0; i < 6; i++)
   memset (D[i], 0, 6 * sizeof (double));

 D[0][0] = 1.0/E1;
 D[1][1] = 1.0/E2; 
 D[2][2] = 1.0/E3; 

 D[0][1] = D[1][0] = -Nu12/E2;
 D[0][2] = D[2][0] = -Nu13/E3;
 D[1][2] = D[2][1] = -Nu23/E3;

 D[3][3] = 1.0/G12;
 D[4][4] = 1.0/G13;
 D[5][5] = 1.0/G23;
}

/* ============================== MathTripBCB ========================= */
/*
** F = B(t)*C*B*const
*/
void MultiplyBCB( int m, int n, double **b, double **c, double cons,
                  double *K )
{
  int    i,j,k,l,iaux;
  double aux1,aux2;
 
  for (i = 0; i < n; i++)
  {
    iaux = i * n;
    for (k = 0; k < m; k++) if (b[k][i] != 0.0)
    {
      aux1 = b[k][i] * cons;
      for (l = 0; l < m; l++) if (c[k][l] != 0.0)
      {
        aux2 = c[k][l] * aux1;
        for (j = i; j < n; j++) K[iaux + j] += b[l][j] * aux2;
      }
    }
  }
}

/* ============================= MathVecMult  ========================= */
void MathVecMult( int n, double **x, double *y, double *z )
{
  int i,j;

  for (i = 0; i < n; i++)
  {
    z[i] = 0.0;
    for (j = 0; j < n; j++) z[i] += x[i][j] * y[j];
  }
}

/* ======================= MathDotVector ============================== */

double MathDotVector(int n, double *a, double *b)
{
  int i;
  double dot = 0.0;

  for (i = 0; i < n; i++) dot += a[i]*b[i];

  return(dot);
}
